## Daniel J. Mallinson and Saahir Shafi
## COVID Policy Response Analysis
## 2020-08-06
## 2020-09-14
## 2020-12-17
## 2021-09-22

rm(list=ls()) #clear workspace

install.packages(c("readxl", "quantmod", "tseries", "plm", "punitroots", "tidyr", "kader"))

library(foreign)
library(readxl)
library(quantmod)
library(tseries)
library(plm) #https://cran.r-project.org/web/packages/plm/vignettes/plmPackage.html
library(reshape2)
library(tidyr)
library(kader)
library(stargazer)

#https://towardsdatascience.com/var-and-panel-data-models-the-powerhouse-of-multivariate-forecasting-techniques-22b8d8888141
#https://www.aptech.com/blog/introduction-to-the-fundamentals-of-panel-data/
#https://cran.r-project.org/web/packages/plm/vignettes/A_plmPackage.html

data <- read.csv("CovidData_9.16.21.csv", fileEncoding="UTF-8-BOM")
data.week <- read.csv("CovidData_2.17.22.csv", fileEncoding="UTF-8-BOM")

data <- data[which(data$Date>20200121),] # Remove observations before Jan 22, 2020

countries <- as.character(unique(data$CountryName))

#add day count to data for plotting
days <- seq(1,345,1)

for(i in 1:length(unique(data$CountryName))){
  temp <- data[which(data$CountryName == countries[i]),]
  temp <- cbind(temp, days)
  names(temp)[length(temp)] <- "day"
  if(i == 1){
    newdata <- temp
  }else{
    newdata <- rbind(newdata, temp)
  }
}

data <- newdata

#Determine countries with missing data
missing <- data[rowSums(is.na(data)) > 0,] #Check for remaining NAs
100-sum(unique(missing$totalpop), na.rm=TRUE)/7800000000*100 #World's population represented by included countries

unique(missing$CountryName)

data <- drop_na(data)
data.week <- drop_na(data.week)

countries <- as.character(unique(data$CountryName))


############################## Figures 1 and 2 ###########################

## Plot Stringency Over Time
#pdf("figure1.pdf", height=8.5, width=11)
tiff("figure1.tiff", height=8.5, width=11, units="in", res=600)
plot.new()
par(fig=c(0,0.5,0,0.5), new=TRUE)
plot.window(xlim=c(0,345), ylim=c(0,100))
for(i in 1:length(countries)){
  country <- countries[i]
  lines(data$day[data$CountryName==country],data$StringencyIndex[data$CountryName==country], col="gray50")
}
lines(data$day[data$CountryName=="United States"],data$StringencyIndex[data$CountryName=="United States"], col="dodgerblue4", lwd=3)
text(150, 20, labels="United States", font=2, col="dodgerblue4")
lines(data$day[data$CountryName=="France"],data$StringencyIndex[data$CountryName=="France"], col="darkred", lwd=3)
text(150, 95, labels="France", font=2, col="darkred")
axis(1, at=c(11, 40, 71, 101, 132, 162, 193, 224, 254, 285, 315), labels=c("2/1", "3/1", "4/1", "5/1", "6/1", "7/1", "8/1", "9/1", "10/1", "11/1", "12/1"))
axis(2, at=seq(0,100,10), labels=seq(0,100,10), las=2)
title(ylab="Stringency Index")
#dev.off()

#pdf("govtresponse_allcountries.pdf", height=4, width=8)
#plot.new()
par(fig=c(0.5,1,0,0.5), new=TRUE)
plot.window(xlim=c(0,345), ylim=c(0,100))
for(i in 1:length(countries)){
  country <- countries[i]
  lines(data$day[data$CountryName==country],data$GovernmentResponseIndex[data$CountryName==country], col="gray50")
}
lines(data$day[data$CountryName=="United States"],data$GovernmentResponseIndex[data$CountryName=="United States"], col="dodgerblue4", lwd=3)
text(150, 20, labels="United States", font=2, col="dodgerblue4")
lines(data$day[data$CountryName=="France"],data$StringencyIndex[data$CountryName=="France"], col="darkred", lwd=3)
text(150, 95, labels="France", font=2, col="darkred")
axis(1, at=c(11, 40, 71, 101, 132, 162, 193, 224, 254, 285, 315), labels=c("2/1", "3/1", "4/1", "5/1", "6/1", "7/1", "8/1", "9/1", "10/1", "11/1", "12/1"))
axis(2, at=seq(0,100,10), labels=seq(0,100,10), las=2)
title(ylab="Government Response Index")
#dev.off()

#pdf("containment_allcountries.pdf", height=4, width=8)
#plot.new()
par(fig=c(0,0.5,0.5,1), new=TRUE)
plot.window(xlim=c(0,345), ylim=c(0,100))
for(i in 1:length(countries)){
  country <- countries[i]
  lines(data$day[data$CountryName==country],data$ContainmentHealthIndex[data$CountryName==country], col="gray50")
}
lines(data$day[data$CountryName=="United States"],data$ContainmentHealthIndex[data$CountryName=="United States"], col="dodgerblue4", lwd=3)
text(150, 20, labels="United States", font=2, col="dodgerblue4")
lines(data$day[data$CountryName=="France"],data$StringencyIndex[data$CountryName=="France"], col="darkred", lwd=3)
text(150, 95, labels="France", font=2, col="darkred")
axis(1, at=c(11, 40, 71, 101, 132, 162, 193, 224, 254, 285, 315), labels=c("2/1", "3/1", "4/1", "5/1", "6/1", "7/1", "8/1", "9/1", "10/1", "11/1", "12/1"))
axis(2, at=seq(0,100,10), labels=seq(0,100,10), las=2)
title(ylab="Containment Health Index")
#dev.off()

#pdf("econsupport_allcountries.pdf", height=4, width=8)
#plot.new()
par(fig=c(0.5,1,0.5,1), new=TRUE)
plot.window(xlim=c(0,345), ylim=c(0,100))
for(i in 1:length(countries)){
  country <- countries[i]
  lines(data$day[data$CountryName==country],data$EconomicSupportIndex[data$CountryName==country], col="gray50")
}
lines(data$day[data$CountryName=="United States"],data$EconomicSupportIndex[data$CountryName=="United States"], col="dodgerblue4", lwd=3)
text(150, 20, labels="United States", font=2, col="dodgerblue4")
lines(data$day[data$CountryName=="France"],data$StringencyIndex[data$CountryName=="France"], col="darkred", lwd=3)
text(150, 95, labels="France", font=2, col="darkred")
axis(1, at=c(11, 40, 71, 101, 132, 162, 193, 224, 254, 285, 315), labels=c("2/1", "3/1", "4/1", "5/1", "6/1", "7/1", "8/1", "9/1", "10/1", "11/1", "12/1"))
axis(2, at=seq(0,100,10), labels=seq(0,100,10), las=2)
title(ylab="Economic Support Index")
dev.off()

#pdf("figure2.pdf", height=8, width=8)
tiff("figure2.tiff", height=8, width=8, units="in", res=600)
plot.new()
par(fig=c(0,1,0,0.5), new=TRUE)
plot.window(xlim=c(0,345), ylim=c(0,18000))
for(i in 1:length(countries)){
  country <- countries[i]
  lines(data$day[data$CountryName==country],data$fiscalmeasures[data$CountryName==country], col="gray50")
}
lines(data$day[data$CountryName=="United States"],data$fiscalmeasures[data$CountryName=="United States"], col="dodgerblue4", lwd=3)
text(150, 10000, labels="United States", font=2, col="dodgerblue4")
lines(data$day[data$CountryName=="France"],data$StringencyIndex[data$CountryName=="France"], col="darkred", lwd=3)
text(150, 2000, labels="France", font=2, col="darkred")
axis(1, at=c(11, 40, 71, 101, 132, 162, 193, 224, 254, 285, 315), labels=c("2/1", "3/1", "4/1", "5/1", "6/1", "7/1", "8/1", "9/1", "10/1", "11/1", "12/1"))
axis(2, at=seq(0,18000,2000), labels=c("0", "2K", "4K", "6K", "8K", "10K", "12K", "14K", "16K", "18K"), las=2)
mtext(2, text="Fiscal Measures ($ Per Capita)", line=3)
#dev.off()

#pdf("healthinvestment_allcountries.pdf", height=4, width=8)
#plot.new()
par(fig=c(0,1,0.5,1), new=TRUE)
plot.window(xlim=c(0,345), ylim=c(0,1300))
for(i in 1:length(countries)){
  country <- countries[i]
  lines(data$day[data$CountryName==country],data$healthinvestment[data$CountryName==country], col="gray50")
}
lines(data$day[data$CountryName=="United States"],data$healthinvestment[data$CountryName=="United States"], col="dodgerblue4", lwd=3)
text(150, 600, labels="United States", font=2, col="dodgerblue4")
lines(data$day[data$CountryName=="France"],data$StringencyIndex[data$CountryName=="France"], col="darkred", lwd=3)
text(150, 200, labels="France", font=2, col="darkred")
axis(1, at=c(11, 40, 71, 101, 132, 162, 193, 224, 254, 285, 315), labels=c("2/1", "3/1", "4/1", "5/1", "6/1", "7/1", "8/1", "9/1", "10/1", "11/1", "12/1"))
axis(2, at=seq(0,1300,200), labels=seq(0,1300,200), las=2)
title(ylab="Health Investment ($ Per Capita)")
dev.off()

########################### Figure 3 #############################

## Response-Risk Ratio
data$totalcaserate <- data$ConfirmedCases/data$totalpop*1000

#pdf("figure3.pdf", height=8.5, width=11)
tiff("figure3.tiff", height=8.5, width=11, units="in", res=600)

cntry <- c("India", "El Salvador", "Trinidad and Tobago", "New Zealand", "United States",
          "Sweden", "Switzerland", "Belarus")

plot.new()
par(fig=c(0,0.5,0,0.5), new=TRUE)
for(i in 1:length(countries)){
  country <- countries[i]
  working <- cbind(data$StringencyIndex[data$CountryName==country], data$totalcaserate[data$CountryName==country])
  max.p <- max(working[,1])
  max.c <- max(working[,2])
  if(i==1){
    max.si <- as.data.frame(cbind(country, max.p, max.c))
  }else{
    max.si <- rbind(max.si, cbind(country, max.p, max.c))
  }
}

max.si[,2] <- as.numeric(paste(max.si[,2]))
max.si[,3] <- as.numeric(paste(max.si[,3]))

max.si <- max.si[order(max.si$max.c),]

lab <- subset(max.si, max.si$country %in% cntry)

plot(max.si$max.p~max.si$max.c, ylab="Maximum Stringency", xlab="Cases Per 1,000", pch=20, col="grey60")
points(lab$max.p~lab$max.c, pch=18, col="black", cex=1.5)
text(lab$max.c, lab$max.p, lab$country, pos=4, cex=.65)
loess.ms <- loess(max.p~max.c, data=max.si)
lines(max.si$max.c, predict(loess.ms), col="blue")

par(fig=c(0,0.5,0.5,1), new=TRUE)

for(i in 1:length(countries)){
  country <- countries[i]
  working <- cbind(data$GovernmentResponseIndex[data$CountryName==country], data$totalcaserate[data$CountryName==country])
  max.p <- max(working[,1])
  max.c <- max(working[,2])
  if(i==1){
    max.gri <- as.data.frame(cbind(country, max.p, max.c))
  }else{
    max.gri <- rbind(max.gri, cbind(country, max.p, max.c))
  }
}

max.gri[,2] <- as.numeric(paste(max.gri[,2]))
max.gri[,3] <- as.numeric(paste(max.gri[,3]))

max.gri <- max.gri[order(max.gri$max.c),]

lab <- subset(max.gri, max.gri$country %in% cntry)

plot(max.gri$max.p~max.gri$max.c, ylab="Maximum Government Response", xlab="Cases Per 1,000", pch=20, ylim=c(0,100), col="grey60")
points(lab$max.p~lab$max.c, pch=18, col="black",cex=1.5)
text(lab$max.c, lab$max.p, lab$country, pos=4, cex=.7)
loess.gri <- loess(max.p~max.c, data=max.gri)
lines(max.gri$max.c, predict(loess.gri), col="blue")

par(fig=c(0.5,1,0,0.5), new=TRUE)

for(i in 1:length(countries)){
  country <- countries[i]
  working <- cbind(data$ContainmentHealthIndex[data$CountryName==country], data$totalcaserate[data$CountryName==country])
  max.p <- max(working[,1])
  max.c <- max(working[,2])
  if(i==1){
    max.chi <- as.data.frame(cbind(country, max.p, max.c))
  }else{
    max.chi <- rbind(max.chi, cbind(country, max.p, max.c))
  }
}

max.chi[,2] <- as.numeric(paste(max.chi[,2]))
max.chi[,3] <- as.numeric(paste(max.chi[,3]))

max.chi <- max.chi[order(max.chi$max.c),]

lab <- subset(max.chi, max.chi$country %in% cntry)

plot(max.chi$max.p~max.chi$max.c, ylab="Maximum Health Containment", xlab="Cases Per 1,000", pch=20, col="grey60")
points(lab$max.p~lab$max.c, pch=18, col="black", cex=1.5)
text(lab$max.c, lab$max.p, lab$country, pos=4, cex=.7)
loess.chi <- loess(max.p~max.c, data=max.chi)
lines(max.chi$max.c, predict(loess.chi), col="blue")

par(fig=c(0.5,1,0.5,1), new=TRUE)

for(i in 1:length(countries)){
  country <- countries[i]
  working <- cbind(data$EconomicSupportIndex[data$CountryName==country], data$totalcaserate[data$CountryName==country])
  max.p <- max(working[,1])
  max.c <- max(working[,2])
  if(i==1){
    max.esi <- as.data.frame(cbind(country, max.p, max.c))
  }else{
    max.esi <- rbind(max.esi, cbind(country, max.p, max.c))
  }
}

max.esi[,2] <- as.numeric(paste(max.esi[,2]))
max.esi[,3] <- as.numeric(paste(max.esi[,3]))

max.esi <- max.esi[order(max.esi$max.c),]

lab <- subset(max.esi, max.esi$country %in% cntry)

cntry1 <- c("New Zealand", "El Salvador", "Belarus")
cntry2 <- c("Trinidad and Tobago", "Sweden", "United States")
cntry3 <- c("Switzerland")
lab1 <- subset(max.esi, max.esi$country %in% cntry1)
lab2 <- subset(max.esi, max.esi$country %in% cntry2)
lab3 <- subset(max.esi, max.esi$country %in% cntry3)

plot(max.esi$max.p~max.esi$max.c, ylab="Maximum Economic Support", xlab="Cases Per 1,000", pch=20, col="grey60")
points(lab$max.p~lab$max.c, pch=18, col="black", cex=1.5)
text(lab1$max.c, lab1$max.p, lab1$country, pos=4, cex=.65)
text(lab2$max.c, lab2$max.p, lab2$country, cex=.65, adj=c(.2,2))
text(lab3$max.c, lab3$max.p, lab3$country, pos=3, cex=.65)
loess.esi <- loess(max.p~max.c, data=max.esi)
lines(max.esi$max.c, predict(loess.esi), col="blue")

dev.off()

############################ Table 1 ###############################

data.week$fiscalmeasures <- log(data.week$fiscalmeasures+1)
data.week$healthinvestment <- log(data.week$healthinvestment+1)
data.week$gdp_norm <- data.week$gdp_norm/1000
data.week$agedpop <- data.week$agedpop*100

# Identify countries with all 0 for economic support, fiscal measures, and health investment
for(i in 1:length(countries)){
  country <- as.character(countries[i])
  working <- data.week$EconomicSupportIndex[data.week$CountryName==country]
  if(sum(working)==0){
    print(country)
  }
}

pdata <- data.week[which(data.week$CountryName!="Belarus"),]
pdata <- pdata[which(data.week$CountryName!="Burkina Faso"),]
pdata <- pdata[which(pdata$CountryName!="Kiribati"),]
pdata <- pdata[which(pdata$CountryName!="Liberia"),]
pdata <- pdata[which(pdata$CountryName!="Libya"),]
pdata <- pdata[which(pdata$CountryName!="Mozambique"),]
pdata <- pdata[which(pdata$CountryName!="Nicaragua"),]
pdata <- pdata[which(pdata$CountryName!="Tanzania"),]

pdata$week <- c(1:46)

## Create pdata.frame for plm
pdata <- pdata.frame(pdata, index=c("CountryName", "week"), drop.index=TRUE, row.names=TRUE)

## Test for unit roots
# Undifferenced
purtest(pdata$casesgrowthfactor, pmax=3, exo="trend", test="hadri", lags="AIC")
purtest(pdata$casesnewweekly, pmax=3, exo="trend", test="hadri", lags="AIC")
purtest(pdata$deathsgrowthfactor, pmax=3, exo="trend", test="hadri", lags="AIC")
purtest(pdata$deathsnewweekly, pmax=3, exo="trend", test="hadri", lags="AIC")
purtest(pdata$StringencyIndex, pmax=3, exo="trend", test="hadri", lags="AIC")
purtest(pdata$GovernmentResponseIndex, pmax=3, exo="trend", test="hadri", lags="AIC")
purtest(pdata$ContainmentHealthIndex, pmax=3, exo="trend", test="hadri", lags="AIC")
purtest(pdata$EconomicSupportIndex, pmax=3, exo="trend", test="hadri", lags="AIC")


# Differenced
pdata$casesgrowthfactor <- diff(pdata$casesgrowthfactor)
pdata$casesnewweekly <- diff(pdata$casesnewweekly)
pdata$deathsgrowthfactor <- diff(pdata$deathsgrowthfactor)
pdata$deathsnewweekly <- diff(pdata$deathsnewweekly)
pdata$StringencyIndex <- diff(pdata$StringencyIndex)
pdata$GovernmentResponseIndex <- diff(pdata$GovernmentResponseIndex)
pdata$ContainmentHealthIndex <- diff(pdata$ContainmentHealthIndex)

pdata <- na.omit(pdata)

purtest(pdata$casesgrowthfactor, pmax=3, exo="trend", test="hadri", lags="AIC")
purtest(pdata$casesnewweekly, pmax=3, exo="trend", test="hadri", lags="AIC")
purtest(pdata$deathsgrowthfactor, pmax=3, exo="trend", test="hadri", lags="AIC")
purtest(pdata$deathsnewweekly, pmax=3, exo="trend", test="hadri", lags="AIC")
purtest(pdata$StringencyIndex, pmax=3, exo="trend", test="hadri", lags="AIC")
purtest(pdata$GovernmentResponseIndex, pmax=3, exo="trend", test="hadri", lags="AIC")
purtest(pdata$ContainmentHealthIndex, pmax=3, exo="trend", test="hadri", lags="AIC")

pdata$StringencyIndex <- diff(pdata$StringencyIndex)
pdata$GovernmentResponseIndex <- diff(pdata$GovernmentResponseIndex)
pdata$ContainmentHealthIndex <- diff(pdata$ContainmentHealthIndex)

pdata <- na.omit(pdata)

purtest(pdata$StringencyIndex, pmax=3, exo="trend", test="hadri", lags="AIC")
purtest(pdata$GovernmentResponseIndex, pmax=3, exo="trend", test="hadri", lags="AIC")
purtest(pdata$ContainmentHealthIndex, pmax=3, exo="trend", test="hadri", lags="AIC")

#### Analysis for Table 1

### Government Response Index

#Fixed Effects models
gri1.fe <- plm(casesgrowthfactor ~  lag(casesgrowthfactor, 1) + lag(GovernmentResponseIndex, 1) + lag(EconomicSupportIndex, 1) + lag(fiscalmeasures, 1) + lag(healthinvestment, 1) + agedpop + gdp_norm + popdensity, data=pdata, model="pooling", effect="twoways")
gri2.fe <- plm(casesnewweekly ~  lag(casesnewweekly, 1) + lag(GovernmentResponseIndex, 1) + lag(EconomicSupportIndex, 1) + lag(fiscalmeasures, 1) + lag(healthinvestment, 1) + agedpop + gdp_norm + popdensity, data=pdata, model="pooling", effect="twoways")
gri1.fe <- plm(deathsgrowthfactor ~  lag(deathsgrowthfactor, 1) + lag(GovernmentResponseIndex, 1) + lag(EconomicSupportIndex, 1) + lag(fiscalmeasures, 1) + lag(healthinvestment, 1) + agedpop + gdp_norm + popdensity, data=pdata, model="pooling", effect="twoways")
gri4.fe <- plm(deathsnewweekly ~  lag(deathsnewweekly, 1) + lag(GovernmentResponseIndex, 1) + lag(EconomicSupportIndex, 1) + lag(fiscalmeasures, 1) + lag(healthinvestment, 1) + agedpop + gdp_norm + popdensity, data=pdata, model="pooling", effect="twoways")

#Random effects test
#H0: Zero variance in individual-specific errors
plmtest(gri1.fe, effect="twoways", type="ghm")
plmtest(gri2.fe, effect="twoways", type="ghm")
plmtest(gri3.fe, effect="twoways", type="ghm")
plmtest(gri4.fe, effect="twoways", type="ghm")

gri1.re <- plm(casesgrowthfactor ~  lag(casesgrowthfactor, 1) + lag(GovernmentResponseIndex, 1) + lag(EconomicSupportIndex, 1) + lag(fiscalmeasures, 3) + lag(healthinvestment, 1) + agedpop + gdp_norm + popdensity, data=pdata, model="random", random.method="walhus")
summary (gri1.re)
gri2.re <- plm(casesnewweekly ~  lag(casesnewweekly, 1) + lag(GovernmentResponseIndex, 1) + lag(EconomicSupportIndex, 1) + lag(fiscalmeasures, 3) + lag(healthinvestment, 1) + agedpop + gdp_norm + popdensity, data=pdata, model="random", random.method="walhus")
summary (gri2.re)
gri3.re <- plm(deathsgrowthfactor ~  lag(deathsgrowthfactor, 1) + lag(GovernmentResponseIndex, 1) + lag(EconomicSupportIndex, 1) + lag(fiscalmeasures, 3) + lag(healthinvestment, 1) + agedpop + gdp_norm + popdensity, data=pdata, model="random", random.method="walhus")
summary (gri3.re)
gri4.re <- plm(deathsnewweekly ~  lag(deathsnewweekly, 1) + lag(GovernmentResponseIndex, 1) + lag(EconomicSupportIndex, 1) + lag(fiscalmeasures, 3) + lag(healthinvestment, 1) + agedpop + gdp_norm + popdensity, data=pdata, model="random", random.method="walhus")
summary (gri4.re)

#Hausmann Test to see if RE improve the model
#HO: Random Effects model consistent and efficient (preferred)
#H1: Fixed effects model consistent and preferred
phtest(gri1.fe, gri1.re) #FE is preferred
phtest(gri2.fe, gri2.re) #FE is preferred
phtest(gri3.fe, gri3.re) #FE is preferred
phtest(gri4.fe, gri4.re) #RE is preferred

## Government Response Index
#pdf("gri.pdf", height=8, width=8)
tiff("gri.tiff", height=8, width=8, units="in", res=600)
#Cases Growth Factor
gri1.output <- list()
for(i in 1:5){
  gri1.fe <- plm(casesgrowthfactor ~  lag(casesgrowthfactor, 1) + lag(GovernmentResponseIndex,i) + lag(EconomicSupportIndex, i) + lag(fiscalmeasures, i) + lag(healthinvestment, i) + agedpop + gdp_norm + popdensity, data=pdata, model="pooling", effect="twoways")
  gri1.output[[i]] <- gri1.fe
  #print(i)
  #print(summary(gri1.fe))
  lag <- i
  coefficient <- coef(summary(gri1.fe))[3,1]
  se <- coef(summary(gri1.fe))[3,2]
  ll <- coef(summary(gri1.fe))[3,1] - 1.96*coef(summary(gri1.fe))[3,2]
  ul <- coef(summary(gri1.fe))[3,1] + 1.96*coef(summary(gri1.fe))[3,2]
  working.data <- data.frame(lag, coefficient, se, ll, ul)
  if(i ==1){
    gri1.plot <- data.frame(lag, coefficient, se, ll, ul)
  }else{
    gri1.plot<- rbind(gri1.plot, data.frame(lag, coefficient, se, ll, ul))
  }
}

plot.new()
par(fig=c(0,0.5,0.5,1), new=TRUE)
plot.window(xlim=c(1,5),ylim=c(-0.2,0.2))
points(gri1.plot$lag, gri1.plot$coefficient, pch=20)
for(i in 1:5){
  lines(c(i,i), c(gri1.plot$ll[i], gri1.plot$ul[i]))
}
abline(h=0)
axis(1, at=seq(0,5,1), labels=seq(0,5,1))
axis(2, at=seq(-0.2,0.2,0.1), labels=seq(-0.2,0.2,0.1), las=2)
title(main="(a) Cases Growth Factor", ylab="Government Response Index", xlab="Lag Length in Weeks")

#New Weekly Cases
gri2.output <- list()
for(i in 1:5){
  gri2.fe <- plm(casesnewweekly ~  lag(casesnewweekly, 1) + lag(GovernmentResponseIndex,i) + lag(EconomicSupportIndex, i) + lag(fiscalmeasures, i) + lag(healthinvestment, i) + agedpop + gdp_norm + popdensity, data=pdata, model="pooling", effect="twoways")
  gri2.output[[i]] <- gri2.fe
  #print(i)
  #print(summary(gri2.fe))
  lag <- i
  coefficient <- coef(summary(gri2.fe))[3,1]
  se <- coef(summary(gri2.fe))[3,2]
  ll <- coef(summary(gri2.fe))[3,1] - 1.96*coef(summary(gri2.fe))[3,2]
  ul <- coef(summary(gri2.fe))[3,1] + 1.96*coef(summary(gri2.fe))[3,2]
  working.data <- data.frame(lag, coefficient, se, ll, ul)
  if(i ==1){
    gri2.plot <- data.frame(lag, coefficient, se, ll, ul)
  }else{
    gri2.plot<- rbind(gri2.plot, data.frame(lag, coefficient, se, ll, ul))
  }
}

par(fig=c(0.5,1,0.5,1), new=TRUE)
plot.new()
plot.window(xlim=c(1,5),ylim=c(-60,90))
points(gri2.plot$lag, gri2.plot$coefficient, pch=20)
for(i in 1:5){
  lines(c(i,i), c(gri2.plot$ll[i], gri2.plot$ul[i]))
}
abline(h=0)
axis(1, at=seq(0,5,1), labels=seq(0,5,1))
axis(2, at=seq(-60,90,10), labels=seq(-60,90,10), las=2)
title(main="(b) New Weekly Cases", ylab="Government Response Index", xlab="Lag Length in Weeks")

#Deaths Growth Factor
gri3.output <- list()
for(i in 1:5){
  gri3.fe <- plm(deathsgrowthfactor ~  lag(deathsgrowthfactor, 1) + lag(GovernmentResponseIndex,i) + lag(EconomicSupportIndex, i) + lag(fiscalmeasures, i) + lag(healthinvestment, i) + agedpop + gdp_norm + popdensity, data=pdata, model="pooling", effect="twoways")
  gri3.output[[i]] <- gri3.fe
  #print(i)
  #print(summary(gri3.fe))
  lag <- i
  coefficient <- coef(summary(gri3.fe))[3,1]
  se <- coef(summary(gri3.fe))[3,2]
  ll <- coef(summary(gri3.fe))[3,1] - 1.96*coef(summary(gri3.fe))[3,2]
  ul <- coef(summary(gri3.fe))[3,1] + 1.96*coef(summary(gri3.fe))[3,2]
  working.data <- data.frame(lag, coefficient, se, ll, ul)
  if(i ==1){
    gri3.plot <- data.frame(lag, coefficient, se, ll, ul)
  }else{
    gri3.plot<- rbind(gri3.plot, data.frame(lag, coefficient, se, ll, ul))
  }
}

par(fig=c(0,0.5,0,0.5), new=TRUE)
plot.new()
plot.window(xlim=c(1,5),ylim=c(-0.1,0.1))
points(gri3.plot$lag, gri3.plot$coefficient, pch=20)
for(i in 1:5){
  lines(c(i,i), c(gri3.plot$ll[i], gri3.plot$ul[i]))
}
abline(h=0)
axis(1, at=seq(1,5,1), labels=seq(1,5,1))
axis(2, at=seq(-0.1,0.1,0.1), labels=seq(-0.1,0.1,0.1), las=2)
title(main="(c) Deaths Growth Factor", ylab="Government Response Index", xlab="Lag Length in Weeks")

#New Weekly Deaths
gri4.output <- list()
for(i in 1:5){
  gri4.fe <- plm(deathsnewweekly ~  lag(deathsnewweekly, 1) + lag(GovernmentResponseIndex,i) + lag(EconomicSupportIndex, i) + lag(fiscalmeasures, i) + lag(healthinvestment, i) + agedpop + gdp_norm + popdensity, data=pdata, model="pooling", effect="twoways")
  gri4.output[[i]] <- gri4.fe
  #print(i)
  #print(summary(gri4.fe))
  lag <- i
  coefficient <- coef(summary(gri4.fe))[3,1]
  se <- coef(summary(gri4.fe))[3,2]
  ll <- coef(summary(gri4.fe))[3,1] - 1.96*coef(summary(gri4.fe))[3,2]
  ul <- coef(summary(gri4.fe))[3,1] + 1.96*coef(summary(gri4.fe))[3,2]
  working.data <- data.frame(lag, coefficient, se, ll, ul)
  if(i ==1){
    gri4.plot <- data.frame(lag, coefficient, se, ll, ul)
  }else{
    gri4.plot<- rbind(gri4.plot, data.frame(lag, coefficient, se, ll, ul))
  }
}

par(fig=c(0.5,1,0,0.5), new=TRUE)
plot.new()
plot.window(xlim=c(0,5),ylim=c(-2,2))
points(gri4.plot$lag, gri4.plot$coefficient, pch=20)
for(i in 1:5){
  lines(c(i,i), c(gri4.plot$ll[i], gri4.plot$ul[i]))
}
abline(h=0)
axis(1, at=seq(0,5,1), labels=seq(0,5,1))
axis(2, at=seq(-2, 2, 1), labels=seq(-2, 2, 1), las=2)
title(main="(d) New Weekly Deaths", ylab="Government Response Index", xlab="Lag Length in Weeks")
dev.off()

## Economic Support Index
#pdf("esi.pdf", height=8, width=8)
tiff("esi.tiff", height=8, width=8, units="in", res=600)
#Cases Growth Factor
for(i in 1:5){
  esi1.fe <- plm(casesgrowthfactor ~  lag(casesgrowthfactor, 1) + lag(GovernmentResponseIndex,i) + lag(EconomicSupportIndex, i) + lag(fiscalmeasures, i) + lag(healthinvestment, i) + agedpop + gdp_norm + popdensity, data=pdata, model="pooling", effect="twoways")
  #print(i)
  #print(summary(esi1.fe))
  lag <- i
  coefficient <- coef(summary(esi1.fe))[4,1]
  se <- coef(summary(esi1.fe))[4,2]
  ll <- coef(summary(esi1.fe))[4,1] - 1.96*coef(summary(esi1.fe))[4,2]
  ul <- coef(summary(esi1.fe))[4,1] + 1.96*coef(summary(esi1.fe))[4,2]
  working.data <- data.frame(lag, coefficient, se, ll, ul)
  if(i ==1){
    esi1.plot <- data.frame(lag, coefficient, se, ll, ul)
  }else{
    esi1.plot<- rbind(esi1.plot, data.frame(lag, coefficient, se, ll, ul))
  }
}

plot.new()
par(fig=c(0,0.5,0.5,1), new=TRUE)
plot.window(xlim=c(1,5),ylim=c(-0.1,0.1))
points(esi1.plot$lag, esi1.plot$coefficient, pch=20)
for(i in 1:5){
  lines(c(i,i), c(esi1.plot$ll[i], esi1.plot$ul[i]))
}
abline(h=0)
axis(1, at=seq(1,5,1), labels=seq(1,5,1))
axis(2, at=seq(-0.1,0.1,0.1), labels=seq(-0.1,0.1,0.1), las=2)
title(main="(a) Cases Growth Factor", ylab="Economic Support Index", xlab="Lag Length in Weeks")

#New Weekly Cases
for(i in 1:5){
  esi2.fe <- plm(casesnewweekly ~  lag(casesnewweekly, 1) + lag(GovernmentResponseIndex,i) + lag(EconomicSupportIndex, i) + lag(fiscalmeasures, i) + lag(healthinvestment, i) + agedpop + gdp_norm + popdensity, data=pdata, model="pooling", effect="twoways")
  #print(i)
  #print(summary(esi2.fe))
  lag <- i
  coefficient <- coef(summary(esi2.fe))[4,1]
  se <- coef(summary(esi2.fe))[4,2]
  ll <- coef(summary(esi2.fe))[4,1] - 1.96*coef(summary(esi2.fe))[4,2]
  ul <- coef(summary(esi2.fe))[4,1] + 1.96*coef(summary(esi2.fe))[4,2]
  working.data <- data.frame(lag, coefficient, se, ll, ul)
  if(i ==1){
    esi2.plot <- data.frame(lag, coefficient, se, ll, ul)
  }else{
    esi2.plot<- rbind(esi2.plot, data.frame(lag, coefficient, se, ll, ul))
  }
}

par(fig=c(0.5,1,0.5,1), new=TRUE)
plot.new()
plot.window(xlim=c(1,5),ylim=c(-30,20))
points(esi2.plot$lag, esi2.plot$coefficient, pch=20)
for(i in 1:5){
  lines(c(i,i), c(esi2.plot$ll[i], esi2.plot$ul[i]))
}
abline(h=0)
axis(1, at=seq(1,5,1), labels=seq(1,5,1))
axis(2, at=seq(-30,20,10), labels=seq(-30,20,10), las=2)
title(main="(b) New Weekly Cases", ylab="Economic Support Index", xlab="Lag Length in Weeks")

#Deaths Growth Factor
for(i in 1:5){
  esi3.fe <- plm(deathsgrowthfactor ~  lag(deathsgrowthfactor, 1) + lag(GovernmentResponseIndex,i) + lag(EconomicSupportIndex, i) + lag(fiscalmeasures, i) + lag(healthinvestment, i) + agedpop + gdp_norm + popdensity, data=pdata, model="pooling", effect="twoways")
  #print(i)
  #print(summary(esi3.fe))
  lag <- i
  coefficient <- coef(summary(esi3.fe))[4,1]
  se <- coef(summary(esi3.fe))[4,2]
  ll <- coef(summary(esi3.fe))[4,1] - 1.96*coef(summary(esi3.fe))[4,2]
  ul <- coef(summary(esi3.fe))[4,1] + 1.96*coef(summary(esi3.fe))[4,2]
  working.data <- data.frame(lag, coefficient, se, ll, ul)
  if(i ==1){
    esi3.plot <- data.frame(lag, coefficient, se, ll, ul)
  }else{
    esi3.plot<- rbind(esi3.plot, data.frame(lag, coefficient, se, ll, ul))
  }
}

par(fig=c(0,0.5,0,0.5), new=TRUE)
plot.new()
plot.window(xlim=c(1,5),ylim=c(-0.03,0.03))
points(esi3.plot$lag, esi3.plot$coefficient, pch=20)
for(i in 1:5){
  lines(c(i,i), c(esi3.plot$ll[i], esi3.plot$ul[i]))
}
abline(h=0)
axis(1, at=seq(1,5,1), labels=seq(1,5,1))
axis(2, at=seq(-0.03,0.03,0.01), labels=seq(-0.03,0.03,0.01), las=2)
title(main="(c) Deaths Growth Factor", ylab="Economic Support Index", xlab="Lag Length in Weeks")

#New Weekly Deaths
for(i in 1:5){
  esi4.fe <- plm(deathsnewweekly ~  lag(deathsnewweekly, 1) + lag(GovernmentResponseIndex,i) + lag(EconomicSupportIndex, i) + lag(fiscalmeasures, i) + lag(healthinvestment, i) + agedpop + gdp_norm + popdensity, data=pdata, model="pooling", effect="twoways")
  #print(i)
  #print(summary(esi4.fe))
  lag <- i
  coefficient <- coef(summary(esi4.fe))[4,1]
  se <- coef(summary(esi4.fe))[4,2]
  ll <- coef(summary(esi4.fe))[4,1] - 1.96*coef(summary(esi4.fe))[4,2]
  ul <- coef(summary(esi4.fe))[4,1] + 1.96*coef(summary(esi4.fe))[4,2]
  working.data <- data.frame(lag, coefficient, se, ll, ul)
  if(i ==1){
    esi4.plot <- data.frame(lag, coefficient, se, ll, ul)
  }else{
    esi4.plot<- rbind(esi4.plot, data.frame(lag, coefficient, se, ll, ul))
  }
}

par(fig=c(0.5,1,0,0.5), new=TRUE)
plot.new()
plot.window(xlim=c(1,5),ylim=c(-0.7,0.1))
points(esi4.plot$lag, esi4.plot$coefficient, pch=20)
for(i in 1:5){
  lines(c(i,i), c(esi4.plot$ll[i], esi4.plot$ul[i]))
}
abline(h=0)
axis(1, at=seq(1,5,1), labels=seq(1,5,1))
axis(2, at=c(-0.7, -0.5, -0.3, -0.1, 0.1), labels=c(-0.7, -0.5, -0.3, -0.1, 0.1), las=2)
title(main="(d) New Weekly Deaths", ylab="Economic Support Index", xlab="Lag Length in Weeks")
dev.off()

## Fiscal Measures
#pdf("fm.pdf", height=8, width=8)
tiff("fm.tiff", height=8, width=8, units="in", res=600)
#Cases Growth Factor
for(i in 1:5){
  fm1.fe <- plm(casesgrowthfactor ~  lag(casesgrowthfactor, 1) + lag(GovernmentResponseIndex,i) + lag(EconomicSupportIndex, i) + lag(fiscalmeasures, i) + lag(healthinvestment, i) + agedpop + gdp_norm + popdensity, data=pdata, model="pooling", effect="twoways")
  #print(i)
  #print(summary(fm1.fe))
  lag <- i
  coefficient <- coef(summary(fm1.fe))[5,1]
  se <- coef(summary(fm1.fe))[5,2]
  ll <- coef(summary(fm1.fe))[5,1] - 1.96*coef(summary(fm1.fe))[5,2]
  ul <- coef(summary(fm1.fe))[5,1] + 1.96*coef(summary(fm1.fe))[5,2]
  working.data <- data.frame(lag, coefficient, se, ll, ul)
  if(i ==1){
    fm1.plot <- data.frame(lag, coefficient, se, ll, ul)
  }else{
    fm1.plot<- rbind(fm1.plot, data.frame(lag, coefficient, se, ll, ul))
  }
}


plot.new()
par(fig=c(0,0.5,0.5,1), new=TRUE)
plot.window(xlim=c(1,5),ylim=c(-0.6,0.6))
points(fm1.plot$lag, fm1.plot$coefficient, pch=20)
for(i in 1:5){
  lines(c(i,i), c(fm1.plot$ll[i], fm1.plot$ul[i]))
}
abline(h=0)
axis(1, at=seq(1,5,1), labels=seq(1,5,1))
axis(2, at=c(-0.6, -0.4, -0.2, 0, 0.2, 0.4, 0.6), labels=c(-0.6, -0.4, -0.2, 0, 0.2, 0.4, 0.6), las=2)
title(main="(a) Cases Growth Factor", ylab="Fiscal Measures (log)", xlab="Lag Length in Weeks")

#New Weekly Cases
for(i in 1:5){
  fm2.fe <- plm(casesnewweekly ~  lag(casesnewweekly, 1) + lag(GovernmentResponseIndex,i) + lag(EconomicSupportIndex, i) + lag(fiscalmeasures, i) + lag(healthinvestment, i) + agedpop + gdp_norm + popdensity, data=pdata, model="pooling", effect="twoways")
  #print(i)
  #print(summary(fm2.fe))
  lag <- i
  coefficient <- coef(summary(fm2.fe))[5,1]
  se <- coef(summary(fm2.fe))[5,2]
  ll <- coef(summary(fm2.fe))[5,1] - 1.96*coef(summary(fm2.fe))[5,2]
  ul <- coef(summary(fm2.fe))[5,1] + 1.96*coef(summary(fm2.fe))[5,2]
  working.data <- data.frame(lag, coefficient, se, ll, ul)
  if(i ==1){
    fm2.plot <- data.frame(lag, coefficient, se, ll, ul)
  }else{
    fm2.plot<- rbind(fm2.plot, data.frame(lag, coefficient, se, ll, ul))
  }
}

par(fig=c(0.5,1,0.5,1), new=TRUE)
plot.new()
plot.window(xlim=c(1,5),ylim=c(-250,250))
points(fm2.plot$lag, fm2.plot$coefficient, pch=20)
for(i in 1:5){
  lines(c(i,i), c(fm2.plot$ll[i], fm2.plot$ul[i]))
}
abline(h=0)
axis(1, at=seq(1,5,1), labels=seq(1,5,1))
axis(2, at=seq(-250,250,50), labels=seq(-250,250,50), las=2)
title(main="(b) New Weekly Cases", ylab="Fiscal Measures (log)", xlab="Lag Length in Weeks")

#Deaths Growth Factor
for(i in 1:5){
  fm3.fe <- plm(deathsgrowthfactor ~  lag(deathsgrowthfactor, 1) + lag(GovernmentResponseIndex,i) + lag(EconomicSupportIndex, i) + lag(fiscalmeasures, i) + lag(healthinvestment, i) + agedpop + gdp_norm + popdensity, data=pdata, model="pooling", effect="twoways")
  #print(i)
  #print(summary(fm3.fe))
  lag <- i
  coefficient <- coef(summary(fm3.fe))[5,1]
  se <- coef(summary(fm3.fe))[5,2]
  ll <- coef(summary(fm3.fe))[5,1] - 1.96*coef(summary(fm3.fe))[5,2]
  ul <- coef(summary(fm3.fe))[5,1] + 1.96*coef(summary(fm3.fe))[5,2]
  working.data <- data.frame(lag, coefficient, se, ll, ul)
  if(i ==1){
    fm3.plot <- data.frame(lag, coefficient, se, ll, ul)
  }else{
    fm3.plot<- rbind(fm3.plot, data.frame(lag, coefficient, se, ll, ul))
  }
}

par(fig=c(0,0.5,0,0.5), new=TRUE)
plot.new()
plot.window(xlim=c(1,5),ylim=c(-.3,.3))
points(fm3.plot$lag, fm3.plot$coefficient, pch=20)
for(i in 1:5){
  lines(c(i,i), c(fm3.plot$ll[i], fm3.plot$ul[i]))
}
abline(h=0)
axis(1, at=seq(1,5,1), labels=seq(1,5,1))
axis(2, at=c(-.3,-.2,-.1,0,.1,.2,.3), labels=c(-.3,-.2,-.1,0,.1,.2,.3), las=2)
title(main="(c) Deaths Growth Factor", ylab="Fiscal Measures (log)", xlab="Lag Length in Weeks")

#New Daily Deaths
for(i in 1:5){
  fm4.fe <- plm(deathsnewweekly ~  lag(deathsnewweekly, 1) + lag(GovernmentResponseIndex,i) + lag(EconomicSupportIndex, i) + lag(fiscalmeasures, i) + lag(healthinvestment, i) + agedpop + gdp_norm + popdensity, data=pdata, model="pooling", effect="twoways")
  #print(i)
  #print(summary(fm4.fe))
  lag <- i
  coefficient <- coef(summary(fm4.fe))[5,1]
  se <- coef(summary(fm4.fe))[5,2]
  ll <- coef(summary(fm4.fe))[5,1] - 1.96*coef(summary(fm4.fe))[5,2]
  ul <- coef(summary(fm4.fe))[5,1] + 1.96*coef(summary(fm4.fe))[5,2]
  working.data <- data.frame(lag, coefficient, se, ll, ul)
  if(i ==1){
    fm4.plot <- data.frame(lag, coefficient, se, ll, ul)
  }else{
    fm4.plot<- rbind(fm4.plot, data.frame(lag, coefficient, se, ll, ul))
  }
}

par(fig=c(0.5,1,0,0.5), new=TRUE)
plot.new()
plot.window(xlim=c(1,5),ylim=c(-6,4))
points(fm4.plot$lag, fm4.plot$coefficient, pch=20)
for(i in 1:5){
  lines(c(i,i), c(fm4.plot$ll[i], fm4.plot$ul[i]))
}
abline(h=0)
axis(1, at=seq(1,5,1), labels=seq(1,5,1))
axis(2, at=seq(-6,4,2), labels=seq(-6,4,2), las=2)
title(main="(d) New Weekly Deaths", ylab="Fiscal Measures (log)", xlab="Lag Length in Weeks")
dev.off()

## Health Investment
#pdf("hi.pdf", height=8, width=8)
tiff("hi.tiff", height=8, width=8, units="in", res=600)
#Cases Growth Factor
for(i in 1:5){
  hi1.fe <- plm(casesgrowthfactor ~  lag(casesgrowthfactor, 1) + lag(GovernmentResponseIndex,i) + lag(EconomicSupportIndex, i) + lag(fiscalmeasures, i) + lag(healthinvestment, i) + agedpop + gdp_norm + popdensity, data=pdata, model="pooling", effect="twoways")
  #print(i)
  #print(summary(hi1.fe))
  lag <- i
  coefficient <- coef(summary(hi1.fe))[6,1]
  se <- coef(summary(hi1.fe))[6,2]
  ll <- coef(summary(hi1.fe))[6,1] - 1.96*coef(summary(hi1.fe))[6,2]
  ul <- coef(summary(hi1.fe))[6,1] + 1.96*coef(summary(hi1.fe))[6,2]
  working.data <- data.frame(lag, coefficient, se, ll, ul)
  if(i ==1){
    hi1.plot <- data.frame(lag, coefficient, se, ll, ul)
  }else{
    hi1.plot<- rbind(hi1.plot, data.frame(lag, coefficient, se, ll, ul))
  }
}

plot.new()
par(fig=c(0,0.5,0.5,1), new=TRUE)
plot.window(xlim=c(1,5),ylim=c(-1,1))
points(hi1.plot$lag, hi1.plot$coefficient, pch=20)
for(i in 1:5){
  lines(c(i,i), c(hi1.plot$ll[i], hi1.plot$ul[i]))
}
abline(h=0)
axis(1, at=seq(1,5,1), labels=seq(1,5,1))
axis(2, at=seq(-1,1,.5), labels=seq(-1,1,.5), las=2)
title(main="(a) Cases Growth Factor", ylab="Health Investments (log)", xlab="Lag Length in Weeks")

#New Weekly Cases
for(i in 1:5){
  hi2.fe <- plm(casesnewweekly ~  lag(casesnewweekly, 1) + lag(GovernmentResponseIndex,i) + lag(EconomicSupportIndex, i) + lag(fiscalmeasures, i) + lag(healthinvestment, i) + agedpop + gdp_norm + popdensity, data=pdata, model="pooling", effect="twoways")
  #print(i)
  #print(summary(hi2.fe))
  lag <- i
  coefficient <- coef(summary(hi2.fe))[6,1]
  se <- coef(summary(hi2.fe))[6,2]
  ll <- coef(summary(hi2.fe))[6,1] - 1.96*coef(summary(hi2.fe))[6,2]
  ul <- coef(summary(hi2.fe))[6,1] + 1.96*coef(summary(hi2.fe))[6,2]
  working.data <- data.frame(lag, coefficient, se, ll, ul)
  if(i ==1){
    hi2.plot <- data.frame(lag, coefficient, se, ll, ul)
  }else{
    hi2.plot<- rbind(hi2.plot, data.frame(lag, coefficient, se, ll, ul))
  }
}

par(fig=c(0.5,1,0.5,1), new=TRUE)
plot.new()
plot.window(xlim=c(1,5),ylim=c(30,550))
points(hi2.plot$lag, hi2.plot$coefficient, pch=20)
for(i in 1:5){
  lines(c(i,i), c(hi2.plot$ll[i], hi2.plot$ul[i]))
}
abline(h=0)
axis(1, at=seq(1,5,1), labels=seq(1,5,1))
axis(2, at=seq(30,550,50), labels=seq(30,550,50), las=2)
title(main="(b) New Weekly Cases", ylab="Health Investments (log)", xlab="Lag Length in Weeks")

#Deaths Growth Factor
for(i in 1:5){
  hi3.fe <- plm(deathsgrowthfactor ~  lag(deathsgrowthfactor, 1) + lag(GovernmentResponseIndex,i) + lag(EconomicSupportIndex, i) + lag(fiscalmeasures, i) + lag(healthinvestment, i) + agedpop + gdp_norm + popdensity, data=pdata, model="pooling", effect="twoways")
  #print(i)
  #print(summary(hi3.fe))
  lag <- i
  coefficient <- coef(summary(hi3.fe))[6,1]
  se <- coef(summary(hi3.fe))[6,2]
  ll <- coef(summary(hi3.fe))[6,1] - 1.96*coef(summary(hi3.fe))[6,2]
  ul <- coef(summary(hi3.fe))[6,1] + 1.96*coef(summary(hi3.fe))[6,2]
  working.data <- data.frame(lag, coefficient, se, ll, ul)
  if(i ==1){
    hi3.plot <- data.frame(lag, coefficient, se, ll, ul)
  }else{
    hi3.plot<- rbind(hi3.plot, data.frame(lag, coefficient, se, ll, ul))
  }
}

par(fig=c(0,0.5,0,0.5), new=TRUE)
plot.new()
plot.window(xlim=c(1,5),ylim=c(-.4,.4))
points(hi3.plot$lag, hi3.plot$coefficient, pch=20)
for(i in 1:5){
  lines(c(i,i), c(hi3.plot$ll[i], hi3.plot$ul[i]))
}
abline(h=0)
axis(1, at=seq(1,5,1), labels=seq(1,5,1))
axis(2, at=seq(-.4,.4,.2), labels=seq(-.4,.4,.2), las=2)
title(main="(c) Deaths Growth Factor", ylab="Health Investments (log)", xlab="Lag Length in Weeks")

#New Weekly Deaths
for(i in 1:5){
  hi4.fe <- plm(deathsnewweekly ~  lag(deathsnewweekly, 1) + lag(GovernmentResponseIndex,i) + lag(EconomicSupportIndex, i) + lag(fiscalmeasures, i) + lag(healthinvestment, i) + agedpop + gdp_norm + popdensity, data=pdata, model="pooling", effect="twoways")
  #print(i)
  #print(summary(hi4.fe))
  lag <- i
  coefficient <- coef(summary(hi4.fe))[6,1]
  se <- coef(summary(hi4.fe))[6,2]
  ll <- coef(summary(hi4.fe))[6,1] - 1.96*coef(summary(hi4.fe))[6,2]
  ul <- coef(summary(hi4.fe))[6,1] + 1.96*coef(summary(hi4.fe))[6,2]
  working.data <- data.frame(lag, coefficient, se, ll, ul)
  if(i ==1){
    hi4.plot <- data.frame(lag, coefficient, se, ll, ul)
  }else{
    hi4.plot<- rbind(hi4.plot, data.frame(lag, coefficient, se, ll, ul))
  }
}

par(fig=c(0.5,1,0,0.5), new=TRUE)
plot.new()
plot.window(xlim=c(1,5),ylim=c(-5,5))
points(hi4.plot$lag, hi4.plot$coefficient, pch=20)
for(i in 1:5){
  lines(c(i,i), c(hi4.plot$ll[i], hi4.plot$ul[i]))
}
abline(h=0)
axis(1, at=seq(1,5,1), labels=seq(1,5,1))
axis(2, at=seq(-5,5,2), labels=seq(-5,5,2), las=2)
title(main="(d) New Weekly Deaths", ylab="Health Investments (log)", xlab="Lag Length in Weeks")
dev.off()

## Containment Health Index
#pdf("chi.pdf", height=8, width=8)
tiff("chi.tiff", height=8, width=8, units="in", res=600)
#Cases Growth Factor
chi1.output <- list()
for(i in 1:5){
  chi1.fe <- plm(casesgrowthfactor ~  lag(casesgrowthfactor, 1) + lag(ContainmentHealthIndex,i) + lag(EconomicSupportIndex, i) + lag(fiscalmeasures, i) + lag(healthinvestment, i) + agedpop + gdp_norm + popdensity, data=pdata, model="pooling", effect="twoways")
  chi1.output[[i]] <- chi1.fe
  #print(i)
  #print(summary(chi1.fe))
  lag <- i
  coefficient <- coef(summary(chi1.fe))[3,1]
  se <- coef(summary(chi1.fe))[3,2]
  ll <- coef(summary(chi1.fe))[3,1] - 1.96*coef(summary(chi1.fe))[3,2]
  ul <- coef(summary(chi1.fe))[3,1] + 1.96*coef(summary(chi1.fe))[3,2]
  working.data <- data.frame(lag, coefficient, se, ll, ul)
  if(i ==1){
    chi1.plot <- data.frame(lag, coefficient, se, ll, ul)
  }else{
    chi1.plot<- rbind(chi1.plot, data.frame(lag, coefficient, se, ll, ul))
  }
}

plot.new()
par(fig=c(0,0.5,0.5,1), new=TRUE)
plot.window(xlim=c(1,5),ylim=c(-.2,.2))
points(chi1.plot$lag, chi1.plot$coefficient, pch=20)
for(i in 1:5){
  lines(c(i,i), c(chi1.plot$ll[i], chi1.plot$ul[i]))
}
abline(h=0)
axis(1, at=seq(1,5,1), labels=seq(1,5,1))
axis(2, at=seq(-.2,.2,.1), labels=seq(-.2,.2,.1), las=2)
title(main="(a) Cases Growth Factor", ylab="Containment Index", xlab="Lag Length in Weeks")

#New Daily Cases
chi2.output <- list()
for(i in 1:5){
  chi2.fe <- plm(casesnewweekly ~  lag(casesnewweekly, 1) + lag(ContainmentHealthIndex,i) + lag(EconomicSupportIndex, i) + lag(fiscalmeasures, i) + lag(healthinvestment, i) + agedpop + gdp_norm + popdensity, data=pdata, model="pooling", effect="twoways")
  chi2.output[[i]] <- chi2.fe
  #print(i)
  #print(summary(chi2.fe))
  lag <- i
  coefficient <- coef(summary(chi2.fe))[3,1]
  se <- coef(summary(chi2.fe))[3,2]
  ll <- coef(summary(chi2.fe))[3,1] - 1.96*coef(summary(chi2.fe))[3,2]
  ul <- coef(summary(chi2.fe))[3,1] + 1.96*coef(summary(chi2.fe))[3,2]
  working.data <- data.frame(lag, coefficient, se, ll, ul)
  if(i ==1){
    chi2.plot <- data.frame(lag, coefficient, se, ll, ul)
  }else{
    chi2.plot<- rbind(chi2.plot, data.frame(lag, coefficient, se, ll, ul))
  }
}

par(fig=c(0.5,1,0.5,1), new=TRUE)
plot.new()
plot.window(xlim=c(1,5),ylim=c(-50,100))
points(chi2.plot$lag, chi2.plot$coefficient, pch=20)
for(i in 1:5){
  lines(c(i,i), c(chi2.plot$ll[i], chi2.plot$ul[i]))
}
abline(h=0)
axis(1, at=seq(1,5,1), labels=seq(1,5,1))
axis(2, at=seq(-50,100,50), labels=seq(-50,100,50), las=2)
title(main="(b) New Weekly Cases", ylab="Containment Health Index", xlab="Lag Length in Weeks")

#Deaths Growth Factor
chi3.output <- list()
for(i in 1:5){
  chi3.fe <- plm(deathsgrowthfactor ~  lag(deathsgrowthfactor, 1) + lag(ContainmentHealthIndex,i) + lag(EconomicSupportIndex, i) + lag(fiscalmeasures, i) + lag(healthinvestment, i) + agedpop + gdp_norm + popdensity, data=pdata, model="pooling", effect="twoways")
  chi3.output[[i]] <- chi3.fe
  #print(i)
  #print(summary(chi3.fe))
  lag <- i
  coefficient <- coef(summary(chi3.fe))[3,1]
  se <- coef(summary(chi3.fe))[3,2]
  ll <- coef(summary(chi3.fe))[3,1] - 1.96*coef(summary(chi3.fe))[3,2]
  ul <- coef(summary(chi3.fe))[3,1] + 1.96*coef(summary(chi3.fe))[3,2]
  working.data <- data.frame(lag, coefficient, se, ll, ul)
  if(i ==1){
    chi3.plot <- data.frame(lag, coefficient, se, ll, ul)
  }else{
    chi3.plot<- rbind(chi3.plot, data.frame(lag, coefficient, se, ll, ul))
  }
}

par(fig=c(0,0.5,0,0.5), new=TRUE)
plot.new()
plot.window(xlim=c(1,5),ylim=c(-0.1,0.1))
points(chi3.plot$lag, chi3.plot$coefficient, pch=20)
for(i in 1:5){
  lines(c(i,i), c(chi3.plot$ll[i], chi3.plot$ul[i]))
}
abline(h=0)
axis(1, at=seq(1,5,1), labels=seq(1,5,1))
axis(2, at=seq(-0.1,0.1,0.1), labels=seq(-0.1,0.1,0.1), las=2)
title(main="(c) Deaths Growth Factor", ylab="Containment Health Index", xlab="Lag Length in Weeks")

#New Daily Deaths
chi4.output <- list()
for(i in 1:5){
  chi4.fe <- plm(deathsnewweekly ~  lag(deathsnewweekly, 1) + lag(ContainmentHealthIndex,i) + lag(EconomicSupportIndex, i) + lag(fiscalmeasures, i) + lag(healthinvestment, i) + agedpop + gdp_norm + popdensity, data=pdata, model="pooling", effect="twoways")
  chi4.output[[i]] <- chi4.fe
  #print(i)
  #print(summary(chi4.fe))
  lag <- i
  coefficient <- coef(summary(chi4.fe))[3,1]
  se <- coef(summary(chi4.fe))[3,2]
  ll <- coef(summary(chi4.fe))[3,1] - 1.96*coef(summary(chi4.fe))[3,2]
  ul <- coef(summary(chi4.fe))[3,1] + 1.96*coef(summary(chi4.fe))[3,2]
  working.data <- data.frame(lag, coefficient, se, ll, ul)
  if(i ==1){
    chi4.plot <- data.frame(lag, coefficient, se, ll, ul)
  }else{
    chi4.plot<- rbind(chi4.plot, data.frame(lag, coefficient, se, ll, ul))
  }
}

par(fig=c(0.5,1,0,0.5), new=TRUE)
plot.new()
plot.window(xlim=c(1,5),ylim=c(-2,2))
points(chi4.plot$lag, chi4.plot$coefficient, pch=20)
for(i in 1:5){
  lines(c(i,i), c(chi4.plot$ll[i], chi4.plot$ul[i]))
}
abline(h=0)
axis(1, at=seq(1,5,1), labels=seq(1,5,1))
axis(2, at=seq(-2,2,1), labels=seq(-2,2,1), las=2)
title(main="(d) New Weekly Deaths", ylab="Containment Health Index", xlab="Lag Length in Weeks")
dev.off()

## Stringency Index
#pdf("si.pdf", height=8, width=8)
tiff("si.tiff", height=8, width=8, units="in", res=600)
#Cases Growth Factor
si1.output <- list()
for(i in 1:5){
  si1.fe <- plm(casesgrowthfactor ~  lag(casesgrowthfactor, 1) + lag(StringencyIndex,i) + lag(EconomicSupportIndex, i) + lag(fiscalmeasures, i) + lag(healthinvestment, i) + agedpop + gdp_norm + popdensity, data=pdata, model="pooling", effect="twoways")
  si1.output[[i]] <- si1.fe
  #print(i)
  #print(summary(si1.fe))
  lag <- i
  coefficient <- coef(summary(si1.fe))[3,1]
  se <- coef(summary(si1.fe))[3,2]
  ll <- coef(summary(si1.fe))[3,1] - 1.96*coef(summary(si1.fe))[3,2]
  ul <- coef(summary(si1.fe))[3,1] + 1.96*coef(summary(si1.fe))[3,2]
  working.data <- data.frame(lag, coefficient, se, ll, ul)
  if(i ==1){
    si1.plot <- data.frame(lag, coefficient, se, ll, ul)
  }else{
    si1.plot<- rbind(si1.plot, data.frame(lag, coefficient, se, ll, ul))
  }
}

plot.new()
par(fig=c(0,0.5,0.5,1), new=TRUE)
plot.window(xlim=c(1,5),ylim=c(-0.1,0.1))
points(si1.plot$lag, si1.plot$coefficient, pch=20)
for(i in 1:20){
  lines(c(i,i), c(si1.plot$ll[i], si1.plot$ul[i]))
}
abline(h=0)
axis(1, at=seq(1,5,1), labels=seq(1,5,1))
axis(2, at=seq(-.1,.1,.1), labels=seq(-.1,.1,.1), las=2)
title(main="(a) Cases Growth Factor", ylab="Stringency Index", xlab="Lag Length in Weeks")

#New Weekly Cases
si2.output <- list()
for(i in 1:5){
  si2.fe <- plm(casesnewweekly ~  lag(casesnewweekly, 1) + lag(StringencyIndex,i) + lag(EconomicSupportIndex, i) + lag(fiscalmeasures, i) + lag(healthinvestment, i) + agedpop + gdp_norm + popdensity, data=pdata, model="pooling", effect="twoways")
  si2.output[[i]] <- si2.fe
  #print(i)
  #print(summary(si2.fe))
  lag <- i
  coefficient <- coef(summary(si2.fe))[3,1]
  se <- coef(summary(si2.fe))[3,2]
  ll <- coef(summary(si2.fe))[3,1] - 1.96*coef(summary(si2.fe))[3,2]
  ul <- coef(summary(si2.fe))[3,1] + 1.96*coef(summary(si2.fe))[3,2]
  working.data <- data.frame(lag, coefficient, se, ll, ul)
  if(i ==1){
    si2.plot <- data.frame(lag, coefficient, se, ll, ul)
  }else{
    si2.plot<- rbind(si2.plot, data.frame(lag, coefficient, se, ll, ul))
  }
}

par(fig=c(0.5,1,0.5,1), new=TRUE)
plot.new()
plot.window(xlim=c(1,5),ylim=c(-40,60))
points(si2.plot$lag, si2.plot$coefficient, pch=20)
for(i in 1:5){
  lines(c(i,i), c(si2.plot$ll[i], si2.plot$ul[i]))
}
abline(h=0)
axis(1, at=seq(1,5,1), labels=seq(1,5,1))
axis(2, at=seq(-40,60,20), labels=seq(-40,60,20), las=2)
title(main="(b) New Weekly Cases", ylab="Stringency Index", xlab="Lag Length in Weeks")

#Deaths Growth Factor
si3.output <- list()
for(i in 1:5){
  si3.fe <- plm(deathsgrowthfactor ~  lag(deathsgrowthfactor, 1) + lag(StringencyIndex,i) + lag(EconomicSupportIndex, i) + lag(fiscalmeasures, i) + lag(healthinvestment, i) + agedpop + gdp_norm + popdensity, data=pdata, model="pooling", effect="twoways")
  si3.output[[i]] <- si3.fe
  #print(i)
  #print(summary(si3.fe))
  lag <- i
  coefficient <- coef(summary(si3.fe))[3,1]
  se <- coef(summary(si3.fe))[3,2]
  ll <- coef(summary(si3.fe))[3,1] - 1.96*coef(summary(si3.fe))[3,2]
  ul <- coef(summary(si3.fe))[3,1] + 1.96*coef(summary(si3.fe))[3,2]
  working.data <- data.frame(lag, coefficient, se, ll, ul)
  if(i ==1){
    si3.plot <- data.frame(lag, coefficient, se, ll, ul)
  }else{
    si3.plot<- rbind(si3.plot, data.frame(lag, coefficient, se, ll, ul))
  }
}

par(fig=c(0,0.5,0,0.5), new=TRUE)
plot.new()
plot.window(xlim=c(1,5),ylim=c(-0.1,0.1))
points(si3.plot$lag, si3.plot$coefficient, pch=20)
for(i in 1:5){
  lines(c(i,i), c(si3.plot$ll[i], si3.plot$ul[i]))
}
abline(h=0)
axis(1, at=seq(1,5,1), labels=seq(1,5,1))
axis(2, at=seq(-0.1,0.1,0.1), labels=seq(-0.1,0.1,0.1), las=2)
title(main="(c) Deaths Growth Factor", ylab="Stringency Index", xlab="Lag Length in Weeks")

#New Daily Deaths
si4.output <- list()
for(i in 1:5){
  si4.fe <- plm(deathsnewweekly ~  lag(deathsnewweekly, 1) + lag(StringencyIndex,i) + lag(EconomicSupportIndex, i) + lag(fiscalmeasures, i) + lag(healthinvestment, i) + agedpop + gdp_norm + popdensity, data=pdata, model="pooling", effect="twoways")
  si4.output[[i]] <- si4.fe
  #print(i)
  #print(summary(si4.fe))
  lag <- i
  coefficient <- coef(summary(si4.fe))[3,1]
  se <- coef(summary(si4.fe))[3,2]
  ll <- coef(summary(si4.fe))[3,1] - 1.96*coef(summary(si4.fe))[3,2]
  ul <- coef(summary(si4.fe))[3,1] + 1.96*coef(summary(si4.fe))[3,2]
  working.data <- data.frame(lag, coefficient, se, ll, ul)
  if(i ==1){
    si4.plot <- data.frame(lag, coefficient, se, ll, ul)
  }else{
    si4.plot<- rbind(si4.plot, data.frame(lag, coefficient, se, ll, ul))
  }
}

par(fig=c(0.5,1,0,0.5), new=TRUE)
plot.new()
plot.window(xlim=c(1,5),ylim=c(-2,1))
points(si4.plot$lag, si4.plot$coefficient, pch=20)
for(i in 1:5){
  lines(c(i,i), c(si4.plot$ll[i], si4.plot$ul[i]))
}
abline(h=0)
axis(1, at=seq(1,5,1), labels=seq(1,5,1))
axis(2, at=seq(-2,1,1), labels=seq(-2,1,1), las=2)
title(main="(d) New Weekly Deaths", ylab="Stringency Index", xlab="Lag Length in Weeks")
dev.off()

#### Supplemental Information

stargazer(gri1.output)
stargazer(gri2.output)
stargazer(gri3.output)
stargazer(gri4.output)
stargazer(chi1.output)
stargazer(chi2.output)
stargazer(chi3.output)
stargazer(chi4.output)
stargazer(si1.output)
stargazer(si2.output)
stargazer(si3.output)
stargazer(si4.output)
